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ABSTRACT 

We simulate the formation and evolution of the local galaxy population starting 
from initial conditions with a smoothed linear density field which matches that derived 
from the IRAS 1.2 Jy galaxy survey. Our simulations track the formation and evolution 
of all dark matter haloes more massive than IO^^Mq out to a distance of 8000 km s^^ 
from the Milky Way. We implement prescriptions similar to those of Kauffmann et al. 
(1999a) to follow the assembly and evolution of the galaxies within these haloes. 
We focus on two variants of the CDM cosmology: a ACDM and a rCDM model. 
Galaxy formation in each is adjusted to reproduce the /-band TuUy-Fisher relation 
of Giovanelli et al. (1997). We compare the present-day luminosity functions, colours, 
morphology and spatial distribution of our simulated galaxies with those of the real 
local population, in particular with the Updated Zwicky Catalog, with the IRAS 
PSCz redshift survey, and with individual local clusters such as Coma, Virgo and 
Perseus. We also use the simulations to study the clustering bias between the dark 
matter and galaxies of differing type. Although some significant discrepancies remain, 
our simulations recover the observed intrinsic properties and the observed spatial 
distribution of local galaxies reasonably well. They can thus be used to calibrate 
methods which use the observed local galaxy population to estimate the cosmic density 
parameter or to draw conclusions about the mechanisms of galaxy formation. To 
facilitate such work, we publically release our z = galaxy catalogues, together with 
the underlying mass distribution. 

Key words: galaxies: clusters: general - galaxies: formation - large-scale structure 
of the Universe 



1 INTRODUCTION 

Over the last decade phenomenological modelling has made 
it possible to follow many aspects of the formation and 
evolution of galaxies within the currently favored hierarchi- 
cal paradigm for the growth of cosmic structure. Recently, 
the grafting of techniques originally developed by White & 
Frenk (1991); Kauffmann et al. (1993); Cole et al. (1994) 
onto high resolution N-body simulations has allowed the spa- 
tial and kinematic distributions of galaxies to be predicted 
in detail as a function of their intrinsic properties (Kauff- 
mann et al. 1997 1999ab; Diaferio et al. 1999 2001; Benson 
et al. 2000ab 2001ab; Springel et al. 2000). This work clar- 
ifies many aspects of the problem of 'galaxy biasing' and 
supersedes the heuristic models previously used to relate 
the galaxy and mass distributions in CDM cosmogonies. Its 
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goal is twofold: to better understand the physical processes 
driving the formation, evolution and clustering of galaxies, 
and to test the standard structure formation paradigm. 

With the steady improvement of computer performance 
and of simulation codes, dissipationless simulations are able 
to achieve ever higher mass resolution. In a recent example, 
Springel et al. (2000, hereafter SOO) followed the formation 
of all galaxies brighter than the Fornax dwarf spheroidal 
within a cluster similar in mass to Coma. In the simulations 
presented below, the use of 7 x 10^ particles allows us to 
follow the formation of all galaxies more luminous than the 
LMC out to 8000 km s"^ from the Milky Way. Improved 
computers have also greatly enhanced the ability of cos- 
mic gas-dynamics codes to simulate galaxy formation (e.g. 
Pearce et al. 2001; Nagamine et al. 2000; White et al. 2000). 
Although such simulations remove the uncertainties due to 
phenomenological modelling of the dynamics and cooling 
of diffuse gas, they retain phenomenological models for the 
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much more uncertain processes of star formation and super- 
nova feedback. Moreover, their greatly increased computa- 
tional cost makes it impossible for them to resolve galaxy 
formation over a volume as large as that modelled in this 
paper. A further major advantage of the techniques used 
here is that the efHciencies for uncertain processes like star 
formation and feedback can be varied to study their influ- 
ence and then adjusted to fit observation without the need 
to run a new simulation for each new parameter set. 

The siunilations we present below use initial conditions 
generated using a technique developed by Kolatt et al. (1996, 
see Bistolas & Hoffman (1998) for a related technique). An 
all-sky redshift survey is smoothed heavily to produce an es- 
timate of the galaxy overdensity field in a spherical volume 
centred on the Milky Way. This is assumed to be a known 
constant times the similarly smoothed local mass overden- 
sity field. One then solves for the linear overdensity field at 
high redshift which would evolve into this quasi-linear local 
field. The initial conditions for the simulation are taken to 
be a random realisation of a Gaussian random field with a 
suitable CDM (or other) power spectrum, but constrained so 
that when suitably smoothed the overdensity field is equal 
to that inferred from the local galaxy distribution. Simula- 
tions run from such initial conditions reproduce the large- 
scale structure of the local universe but have the character- 
istics of the assumed CDM model on small scales where the 
smoothed local galaxy distribution imposes no constraints. 

Although galaxy formation should occur in such simu- 
lations exactly as in random realisations of the underlying 
CDM model, there are a number of advantages to having the 
large-scale structure of the simulation correspond in detail 
to that of the local universe. For certain types of galaxies, in 
particular dwarfs, surveys are restricted to our local neigh- 
borhood. For many types of galaxies, surveys are most com- 
plete and the properties of the galaxies best characterised 
in this region. When interpreting such surveys one must bo 
wary of biases introduced by the particular structure of our 
neighborhood. Such biases are clearly minimised in mod- 
els which reproduce the local structure. Distances to galax- 
ies can only be measured with sufficient accuracy to esti- 
mate their peculiar motions out to redshifts of about 10,000 
km s^^. As a result, detailed studies of large-scale flows are 
only possible in our local neighborhood. Such studies aim 
to verify that flows axe gravitationally induced and to use 
them to measure the cosmic density parameter. Simulations 
of the nearby universe are ideal for calibrating such studies 
and for checking that they produce unbiased estimates of 

The present paper is organized as follows. In Section 
2 we describe both how we construct constrained initial 
conditions from the density field of the IRAS 1.2 Jy sur- 
vey and how we carry out daxk matter simulations from 
these initial conditions. Section 3 explains how halo cata- 
logues and halo merging trees arc built from the simulation 
outputs and summarises the phcriomenological treatment of 
galaxy formation which we graft onto these trees, emphasis- 
ing points where it differs from the treatment in Kauffmann 
et al. (1999a, hereafter K99). This section also compares our 
dark matter distributions to those found in unconstrained 
simulations and our simulated distributions of galaxy lumi- 
nosity, colour and morphology to those observed in large 
surveys. Section 4 begins our detailed comparison with the 



local universe by matching simulated rich clusters object 
by object with real rich clusters. In Section 5 we explain 
how we generate 'mock' catalogues for direct comparison 
with the IRAS PSCz and UZC surveys. In Section 6 we 
use these catalogues to carry out a point-by-point compar- 
ison of the smoothed galaxy and mass density fields, while 
Section 7 extends the comparison to smaller scales using 
cross-correlation statistics. Section 8 displays mock versions 
of the Mark III catalogue of galaxy peculiar velocities to 
illustrate the application of our simulations to cosmic flow 
problems. Finally, we give a brief summary and discussion 
of our results in Section 9. 



2 DARK MATTER 

Our modelling is based on large high resolution simulations 
of the evolution of the dark matter distribution in a spherical 
region surrounding the Milky Way. The first half of this 
section explains how we obtain initial conditions for these 
simulations; the second half describes the code used to follow 
their evolution. 

2.1 Constructing initial conditions 

We use the techniques developed by Kolatt et al. (1996, 
hereafter K96) to set up initial conditions such that at z = 
the evolved mass overdensity in the simulation is a suitably 
scaled version of the galaxy overdensity in the IRAS 1.2 Jy 
survey of Fisher et al. (1994 1995) once both are smoothed 
on the same large scale. We take all the observed galaxies out 
to a redshift of 12,000 km s"^, weight each by the inverse of 
the survey selection function at its redshift, and smooth with 
a Gaussian of 1-D dispersion 5h~^ Mpc to obtain the galaxy 
density field at all points within a sphere of radius 10,000 
km s^^. When smoothing wc take care to account properly 
for the regions of space not sampled by the observational 
survey. These are primarily behind the Galactic Plane and 
in the strip which was not scanned by the IRAS satellite. 

For each cosmology (ACDM and tCDM as in K99) 
we obtain a 'target' mass overdensity field by scaling the 
galaxy overdensities so that their rms within the region out 
to 10,000 km s"'^ is equal to the value obtained by smooth- 
ing the z = linear power spectrum of the particular cos- 
mological model with a Gaussian of dispersion 5 Mpc. 
As described by K96 the Eulerian Zel'dovich-BernouUi equa- 
tion can then be integrated back in time to give the linear 
density field at a = 50 which gives rise to the target smooth 
density field at the present day. This linear field is then 
"Gaussianised" , i.e. it is mapped onto a new field such that 
the new density at each point is a monotonic function of 
the original overdensity and is Gaussian on the observed re- 
gion with rms equal to that expected at « = 50 from the 
smoothed theoretical power spectrum. 

We now apply the Hoffman- Riback algorithm (Hoffman 
& Ribak 1991 1992; Canon & Hoffman 1993) to generate an 
initial displacement field on a 256^ grid which (i) IS peri- 
odic on a simulation cube of side L = 240 Mpc, (ii) 
is a Gaussian random field with power per mode equal to 
that expected at a = 50 for the chosen cosmology for all 
wavenumbers between the fundamental, ko = 2it/L and 
64 ko, and (iii) is constrained so that when smoothed with 
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a Gaussian of comoving dispersion 5 h^^ Mpc it reproduces 
the target initial density field everywhere within a sphere of 
radius 80 Mpc centred on the cube. 

We supplement this constrained low-frequency displace- 
ment field with an unconstrained high-frequency displace- 
ment field constructed as follows. On a periodic cube of 
side 80 Mpc we set up an unconstrained Gaussian ran- 
dom field with power per mode equal to that expected in 
the chosen cosmology for all wavenumbers between 64 feo 
and 2'K/{0.7h^^ Mpc) = 343 ko- The latter corresponds to a 
wavelength equal to twice the mean interparticle separation 
of our high resolution region (see below). We replicate this 
smaller cube 27 times to obtain a high-frequency displace- 
ment field everywhere within the 240 h^^ Mpc cube. 

Because the version of the integration code we use 
has vacuum boundary conditions, it works best on near- 
spherical regions. We therefore create an unperturbed but 
variable resolution particle load within a sphere of radius 
240/1"^ Mpc X = 207. 85 h-^ Mpc. We first fill the 240 

Mpc cube with a uniform but irregular "glass" of equal 
mass particles. (See White (1995) for a discussion of such 
glass distributions.) We then border this cube with periodic 
replications and select the sphere which just encloses the 
central cube as the final region to be simulated. Within this 
large sphere, we excise all particles which will end up less 
than 80 Mpc from the center by z = 0, and we replace 
them with a region of the same size and shape cut from 
a glass distribution with a mean interparticle separation of 
0.35 Mpc. The particle masses are adjusted to ensure 
the correct mean mass density in each of the two regions. 
Table 1 lists the cosmologicaJ parameters and the particle 
numbers and masses for each of our two simulations. 

The final initial condition (at z = 50) is created by in- 
terpolating the high- and low-frequency displacement fields 
to the unperturbed position of each particle, displacing it, 
and assigning it a peculiar velocity proportional to its dis- 
placement in accordance with usual Zcl'dovich approxima- 
tion. We move the low-mass, high-resolution particles in the 
inner region using the sum of the two displacement fields. 
The high-mass, low-resolution particles in the outer region 
are perturbed using the low-frequency constrained displace- 
ment field only. The latter field must be periodically repli- 
cated in order to displace the low resolution particles outside 
the 240 Mpc central cube. 

2.2 Simulating the Dark Matter 

We carried out simulations from the above initial condi- 
tions using the parallel tree-code GADGET (Springel et al. 
2001). Both simulations were started with fixed comoving 
softenings of 70 kpc and 420 kpc for high- and 
low-resolution particles respectively. Once the correspond- 
ing physical comoving softenings reached 20 kpc and 
120 kpc they were kept constant at these values. 

For the ACDM simulation, this resulted in particles in 
the densest regions of the high resolution region requiring 
roughly 75000 adaptive timesteps to reach redshift zero. The 
number of equivalent full timesteps (the total number of 
forces computed divided by the number of particles) was 
much smaller however, of order 400. Both runs were per- 
formed on 512 processors of the CRAY T3E at the Com- 
puter Center (RZG) of the Max-Planck Society at Garching. 
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Figure 2. Mass function of haloes in the simulations. The solid 
and dotted lines are the simulated ACDM and rCDM models 
respectively. The dashed and dashed-dotted lines are the ACDM 
and rCDM predictions of Jenkins et al. (2001). 



The ACDM simulation required about 104 processor-hours 
to complete (the total CPU time consumption was about 
50000 hours), while the rCDM simulation was somewhat 
faster. We saved the particle positions and velocities at 41 
times logarithmically spaced in expansion parameter from 
z = 10.3 to 2 = 0. 

In Fig. 1 we show the present-day dark matter distri- 
bution in our ACDM and rCDM models. The projected re- 
gion is a 30 h^^ Mpc thick slice parallel to the supergalactic 
plane and extending from -15 Mpc to 15 Mpc in 
SGZ. The side of the sUce is 180 Mpc long, and the 
Milky- Way is at its centre. Well-known local structures are 
easily identified, in particular, the Great Attractor region, 
the Coma-Great Wall structure, the Pisces-Perseus super- 
cluster filament, and the Local Void. 

Many of the known nearby rich clusters, including 
Virgo, can also be identified, even though the smoothing 
kernel used to generate the initial conditions is substantially 
larger than a rich cluster and such clusters axe actually quite 
poorly represented in the IRAS 1.2 Jy survey. 

One can start to see the low resolution region at the 
corners of these images which extend more than 80 h^^ Mpc 
from the centre. In the remainder of this paper we shall 
concentrate on analysis of the region within 80 Mpc 
and we will exclude from consideration the few dark haloes 
near the boundary which are contaminated with high-mass 
particles. 

As a test that our simulations are correctly reproducing 
the small scale structure expected in each cosmology, Fig. 2 
compares their halo mass functions at z = with those 
expected according to the large ensemble of unconstrained 
simulations analysed by Jenkins et al. (2001). As discussed 
in the next section, we identify haloes in our simulations 
using a friends-of-friends group-finder with 6 = 0.2 (Davis 
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Table 1. Simulation parameters of our dissipationless cosmological simulations. Nhr, Ni,., Mhr, and rgoft.hn ''soft.lr 
give the number, the mass (units Mq) and the physical softening lengh (units kpc) of the high-resolution and 

low-resolution particles. rosj)Octi\'oly. 



Model 


hioo 


Ho 


A 








^soft,hr 




Ml, 




ACDAl 


0.7 


(),:-! 


0,7 


0.!) 


51 X 10'' 


o.:-!() X 10^" 


20 


20.5 X 10'' 


1. 1 X 10^^ 


120 


rCDM 


0.5 


1.0 


0.0 


0.6 


53 X 10« 


1.2 X 10^° 


20 


20.4 X 10^ 


4.8 X 10" 


120 



Figure 1. The present-day distribution of the dark matter. ACDM and rCDM arc shown on left and right, respectively. The slice shown 
is 180 Mpc wide, 30 Mpc thick and encompasses the supergalactic plane. Prominent clusters are labelled. A high resolution copy 
of this Figure can be found at http://www.mpa-garching.mpg.de/NumCos/CR/High-res/index.html 



et al. 1985). The analytic fit which .Jenkins ot al. (2001) 
use to summarise their data for such FOF mass functions 
is also an excellent description of our own halo data. Note 
that no parameters were adjusted in making this compari- 
son. Significant discrepancies occur only for the largest mass 
haloes. Those arc plausibly a consequence of the constraints 
we impose on our initial conditions. 



3 PUTTING IN THE GALAXIES 

We follow the assembly, the evolution and the merging 
of galaxies within our dark matter simulations using phe- 
nomenological prescriptions based closely on those of K99. 
In this section we give an outline description of the scheme 
as a whole. We then discuss particular aspects of the mod- 
elling where we have found it advantageous to modify the 
K99 prescriptions. Finally we describe how we set the pa- 
rameters which describe uncertain physical processes (for 
example, the efficiencies of star formation and feedback) us- 
ing observed properties of galaxies such as the TuUy-Fisher 
relation and the luminosity function. 

Our galaxy formation scheme is based on the idea, orig- 
inally introduced by White & Rees (1978), that galaxies re- 
sult from the condensation of diffuse gas at the centres of 
dark haloes. As a result all galaxies arc found within dark 
haloes, either as the principal central object or as satellites 
in orbit about it. The first step is thus to identify all the dark 
haloes in each snapshot of the simulation and to construct a 
tree which describes how they merge together from snapshot 
to snapshot. The galaxy population can be reconstructed by 
following evolution down the tree from the earliest times to 
the present. Within each halo at each time diffuse gas cools 
onto the disk of the central galaxy, stars form from cold gas 
both in this galaxy and in its satellites, feedback from this 
star formation reheats some of the cold gas, and satellites 
occasionally merge with the central galaxy, modifying its 
morphology if they are sufficiently massive. Galaxy motions 
arc followed in a dynamically consistent way by attaching 
each galaxy to a dark matter particle. For a central galaxy 
this identification is with the most bound particle of the 
halo. The particle identification remains fixed if this halo 
merges with a more massive one so that the galaxy becomes 
a satellite. 

We now briefly discuss aspects of this modelling, em- 
phasising aspects where our implementation differs in detail 
from that of K99. 



3.1 Tree structure 

As in K99 we build a halo catalogue for each of the 41 snap- 
shots of each simulation using a fricnds-of-friends algorithm 
with linking length 0.2 times the mean interparticle separa- 
tion (Davis et al. 1985). For further analysis we consider only 
haloes with 10 or more particles; smaller haloes frequently 
disappear in later snapshots. The more massive haloes are 
linked from each snapshot to the next to form a tree exactly 
as in K99. Note that although we can pinpoint the location 
and mass of haloes quite well down to our limit of 10 parti- 
cles, we can only reliably track the merging history of haloes 
more massive than about 100 particles. Lower mass haloes 
can typically be traced back through a few snapshots only 
before they vanish below the resolution limit. As a result, 
we can get reasonable estimates for the mass of the cen- 
tral galaxy in a 10 to 20 particle halo, but we cannot infer 
anything about its morphology for halo masses below 100. 

This leads us to define the luminosity resolution limit 
of our simulated galaxy catalogues to be the typical lumi- 
nosity of the central galaxy in a 10 particle halo, and the 
morphology resolution limit to be the typical luminosity of 
the central galaxy in a 100 particle halo. At « = the cir- 
cular velocities of 10 and 100 particle haloes are 45 km 
and 103 km s~^ respectively in our ACDM simulation. For 
rCDM the corresponding numbers arc 76 km s~^ and 160 
km s~^. (We define the circular velocity V200 of a halo to be 
the value measured at R20Q , the radius of the sphere centred 
on the most bound particle within which the mean density is 
200 times the critical density of the Universe.) In the galaxy 
catalogues we construct below, the corresponding luminos- 
ity and morphology limits are Mb = —16.27 and —18.46 for 
ACDM, Mb = -18.50 and -20.60 for tCDM. In this paper, 
"resolution limit" without a qualifier will always denote the 
luminosity resolution limit. 

Galaxies arc of three types in our evolution scheme: cen- 
tral galaxies which are always identified with the most bound 
particle (MBP) of their halo, satellite galaxies which are 
identified with other halo particles, and field galaxies whose 
particle is not a member of any catalogued halo. When set- 
ting up the galaxy population at redshift zo, we adopt the 
following procedures which differ slightly from those of K99. 

(1) If no particle in a halo was a galaxy in the previous 
snapshot at zi > zo, then we create a new galaxy with de- 
fault properties based on a simple infall model for the earlier 
evolution of the halo up to zo- 
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(2) If no particle in the halo was the central galaxy of a 
halo at zi, we identify the most massive (in terms of total 
stellar and cold gas mass) of the galaxies which do populate 
the halo as its central galaxy, we transfer the properties of 
this galaxy to the halo MBP, and wc delete the particle 
previously assigned to it from the galaxy list. 

(3) If at least one particle was a central galaxy at zi, 
we identify the new central galaxy as the one whose halo 
donated the most mass to the zo halo. Its properties are 
transferred to the MBP and its previous particle identifica- 
tion is deleted from the gala^xy list. The other 'old' central 
galaxies become satellite galaxies in the new halo. 

(4) Galaxies from zi which arc part of no halo at zo arc 
designated as field galaxies. These objects are always few in 
number and are near the resolution limit of the simulation. 
They play no significant role in the results we present below. 

This scheme is designed to ensure that a halo central 
galaxy always inherits the properties of the most plausible 
(usually the most massive) candidate for its progenitor in 
the previous snapshot. We have tested a variant in which 
inheritance is based on luminosity rather than mass and find 
it to make no discernable difference. The scheme corrects 
a loophole in the K99 procedures which allowed relatively 
massive haloes occasionally to find no progenitor for their 
central galaxies. As a result they were assigned relatively 
low-mass central galaxies and they appeared as outliers in 
the Tully-Fisher relation. 



3.2 Cooling, star formation and feedback 

The 41 snapshots stored from z = 10.42 to z = sample the 
evolution of the system more coarsely than is desirable when 
simulating the evolution of the galaxy population. As in K99 
we circumvent this by dividing the time between each pair 
of snapshots into 50 small steps, which wo use for integrat- 
ing the simple phenomenological equations which describe 
gas infall and cooling, star formation, feedback and satel- 
lite merging. We have checked that doubling the number of 
steps does not change our results. 

We treat cooling using exactly the same prescriptions 
as K99. These assume that gas cools out of the hot atmo- 
sphere of a halo onto the disk at its centre at a rate which 
depends only on redshift, on halo mass and on halo hot gas 
content. For haloes in which cooling is rapid (at high red- 
shift, relatively small mass, and large hot gas fraction) the 
atmosphere is assumed to condense out on the dynamical 
time of the halo. When cooling times are longer, the con- 
densation rate is taken from a simple cooling fiow model. 
We calculate cooling times using the coUisional ionisation 
cooling curve of Sutherland & Dopita (1993) assuming pri- 
mordial abundances. We do not attempt to follow chemical 
enrichment. We treat the global baryon fraction of the mod- 
els /bar = as an adjustable parameter. 

These prescriptions give rise to a problem that has been 
noticed at least since the work of White & Frenk (1991). In 
galaxy cluster haloes they predict strong cooling flows (sim- 
ilar to those apparently observed in clusters like Perseus) 
which deposit large amounts of cold gas at cluster centres. 
The same star formation prescriptions which produce rea- 
sonable galaxies at the centre of IO^^Mq haloes then predict 
massive, luminous, star-forming galaxies which are quite in- 



consistent with the observed properties of the central galax- 
ies in clusters. This is just another facet of the well-known 
problem of the fate of the gas apparently deposited by cool- 
ing flows (Fabian et al. 1991; AUen & Fabian 1997). 

A number of authors have tried to deal with this prob- 
lem in the present context by building more elaborate mod- 
els for the structure and thermodynamic history of the gas in 
dark matter haloes (Somerville & Primack 1999; Cole et al. 
2000, hereafter SP99, COO). We have tried implementing the 
prescriptions advocated by SP99 and flnd that in our own 
models they produce only a rather limited improvement. As 
a result, we have preferred to retain the simpler but cruder 
ad hoc solution of K99 - wc switch off all cooling of gas in 
haloes with circular velocities exceeding 350 km s^^. This is 
equivalent to assuming that gas processed through cluster- 
like cooling flows does not end up being available for normal 
star formation. It also ensures that no disks are formed with 
circular velocities exceeding the largest observed for real disk 
galaxies. 

Our models for star formation and feedback are also 
identical to those of K99. The star formation rate in a galaxy 
is assumed to be proportional to its cold gas mass divided 
by its dynamical time: 

M. = = 100 a H{z) Mcoid (1) 

^■dyn 

where a is a free parameter to be fitted against observation, 
H{z) is the Hubble parameter, and the second equality fol- 
lows from a very simple model for the typical dynamical 
time of galaxy disks, taken to be 0.1 -R2()()/^^2()0. We use the 
population synthesis models of Bruzual & Chariot (1993) 
to compute the luminosities of our galaxies from their star 
formation histories assuming a Scalo IMF. 

Feedback is assumed to act by reheating some of the 
cold gas in the disk. We write the reheating rate as: 



A/reheat = £ 



4 MtJTsNEsj. 

^ 200 



(2) 



where r/sN = 5 X 10~^Mg^ is the number of supernovae 
produced per solar mass of stars in a Scalo IMF, Esn ~ 
10^^ erg is the mean energy released per supernova, V200 is 
the circular velocity of the halo at R200, and e is a free 
parameter to be tuned to match observations. Since it is 
unclear whether the reheated gas will remain in the halo 
or will be expelled, a model of each type was discussed by 
K99. We have tried both schemes in our own simulations 
and, as in K99, we find that ACDM fits the observations 
better with retention feedback. For rCDM , however, we 
are able to fit the Tully-Fisher relation with less scatter 
using the ejection scheme. This may appear surprising in 
view of fig. 6 of K99, but it is accounted for by our different 
prescription for computing Vdiak (see below). 

3.3 Dust 

We use the simple model of dust extinction proposed by 
SP99. We compute the face-on _B-band optical depth as: 



tb = tb,* 



(3) 



with Ts,« = 0.8, Lb,* = 6 x IO^Lq and (3 = 0.5, as given by 
Wang & Heckman (1996). We then use the galactic extinc- 
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tion curve of Cardelli et al. (1989) to calculate the optical 
depth in other bands. We pick a random inclination for each 
galaxy, and apply a dust correction to its disk component 
only. This is calculated using a "slab" model as: 

(1 _ p-TX sec e \ 
^ (4) 
rAsec6i J 

K99 apply a dust correction only to the disks of galaxies with 
star formation rates SFR > 0.5MQ/yr. Here, we correct all 
disks, regardless of their SFR. In practice, the differences 
induced by the change in prescription are not large, and 
although the new prescription gives somewhat better results, 
it is still unable to eliminate the apparent excess of galaxies 
at the bright end of the luminosity function. 

3.4 Merging 

As in K99, wc allow a satellite galaxy to merge only with 
the central galaxies of its halo. Each satellite galaxy carries a 
merging time counter which is initialized with an estimate of 
the timoscalc for orbital decay through dynamical friction. 
This counter is reinitialized each time the satellite's halo 
undergoes a major merger or is accreted by a more mas- 
sive halo. Otherwise it decreases steadily and the satellite is 
assumed to merge with the central galaxy when it reaches 
zero. 

The initial timescale is based both on analytic estimates 
(Binney & Tremaine 1987) and on fits to numerical simular 
tions (Navarro et al. 1995). Some problems with these for- 
mulae have been pointed out by 800, who note that when 
haloes of comparable mass merge the inferred merging times 
for their central gala^xies are typically much shorter than 
the time actually taken in high resolution simulations of 
the process. This could, in principle, lead to an ovcrmcrging 
problem. There is also disagreement in the literature on the 
appropriate formula for the Coulomb logarithm when esti- 
mating dynamical friction timescales. Somerville & Primack 
(1999) take In A = ln(l-F(Mhaio/Msat)^) while K99 and SOD 
use InA = ln(Mhaio/Msat) and InA = In(l + Mhaio/Msat) 
respectively. Here, we adopt the formula of SP99 since it 
improves the bright-end behaviour of our luminosity func- 
tions. 

We attempt to avoid any overmerging problem by en- 
suring that galaxies never merge on timescales shorter than 
their halo crossing times. If tbic is the merging timescale 
derived from dynamical friction considerations, we initialise 
the merging time counter using: 

^merging — tdyn(2 -\- - '■ ) (5) 
^dyn 

This encapsulates the physical assumption that the satel- 
lite cannot merge with the central galaxy fewer than two 
crossing times after it first enters the halo. 

When the merging time counter vanishes, we add the 
properties of the satellite to those of the central galaxy, and 
we remove the satellite from the galaxy list. We introduce 
a bulge formation threshold /bulge as a free parameter, used 
to match the observed abundance of galaxies by morphol- 
ogy. Mergers in which Msat/Afccntrai > /bulge are considered 
major mergers and create an elliptical gala^xy containing 
the stars of both progenitors, as well as stars formed from 
their cold gas component in a burst assumed to last 10* 



years. A merger with Msat/Mcontrai < /bulge is considered 
minor merger. We add the cold gas of the satellite to that of 
the central galaxy. If no bulge already exists in the central 
galaxy, we create a new bulge with the stars of the merging 
satellite, otherwise we simply add the stars to the existing 
bulge. Such a scheme is sufficient for our current purposes, 
more realistic treatments include that of Somerville et al. 
(2001). 

3.5 Setting parameters 

The models set out above have four free parameters: the 
baryon fraction /bar, the star formation efficiency a, the 
supernova feedback efficiency e, and the bulge formation 
threshold /bulge- We choose these parameters within the 
range which seems physically plausible in such a way as to 
maximise the agreement of the resulting galaxy populations 
with observation. Notice that this procedure is only feasible 
in a scheme of the type we have implemented where changes 
in such physical parameters do not require the supercom- 
puter simulations to be repeated. We attempt to match the 
following observations: 

(1) The velocity-luminosity relation for simulated spirals 
before dust correction should match the published /-band 
relation of Giovanelli et al. (1997), which is already corrected 
for internal extinction. 

(2) The bj-hand luminosity function should be in rough 
agreement with those of the 2dF and SDSS galaxy surveys 
(Folkes et al. 1999; Blanton et al. 2001). Note that our poor 
resolution prevents meaningful comparison of the faint end 
slopes so we concentrate on getting reasonable fits for the 
brighter galaxies. 

(3) Galaxies with Vdisk ~ 220 km s~^ should host ~ 
IO^^Mq of stars and a few times lO^M© of cold gas (c.f. 
K99, SP99). The latter includes both molecular and atomic 
gas. 

(4) The morphology distribution should agree with the 
observation: this is discussed further in section 3.5.4. 

In Table 2 we list our estimated best parameters given 
these constraints, which are described in more detail in later 
subsections. Note that the baryon fraction in the tCDM 
model is high and inconsistent with standard big-bang 
nucleosynthesis (Buries et al. 1999), although it matches 
the baryon fraction observed in massive clusters (Ettori & 
Fabian 1999). Interestingly, if we adopt the preferred cool- 
ing model of SP99, we can fit the Tully-Fisher relation with 
baryon fractions of order /bar 0.08 and /bar 0.15 for 
ACDM and rCDM respectively, albeit with overbright "cD" 
galaxies. This same trend in the required /bar was also found 
by SP99. Note that taking into account the chemical enrich- 
ment of the hot halo gas would boost the cooling rate at 
late times, and allow the observed Tully-Fisher relation to 
be matched for a lower baryon fraction. (See, for example, 
models n and c of van Kampen et al. 1999.) We also find 
that lower baryon fractions are allowed if we assume a sub- 
stantial fraction (~ 20%) of the mass of newly formed stars 
is returned to the cold gas of the galaxy, as expected from 
SNe and stellar winds (sec COO). Wc did not include this 
possibility in the models described here, as it makes more 
sense to use it in combination with a full chemical enrich- 
ment scheme. 
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Figure 3. The z=0 galaxy distribution in ACDM (left panel) and rCDM (right panel) cosmologies. The region shown is the same as in 
Fig. 1. We plot all galaxies brighter than Mb < —19.5 and Mb < —21.3 in ACDM and tCDM respectively. This criterion selects ~ 3500 
galaxies in both cases. The size of the symbols scales with the B-band luminosity and the colours follow B — V index. A high resolution 
copy of this Figure can be found at http://www.mpa-garching.mpg.de/NumCos/CR/High-res/index.htmI 



Table 2. Parameters of the galaxy formation models, /b^r = 
^bar/f^o is the baryon fraction, /bulge the bulge formation thresh- 
old, a and e the star formation and feedback efficiencies respec- 
tively, and "feedback" gives the scheme that we have used. 



Model 


/bar 


a 


e 


/bulge 


feedback 


ACDM 


0.12 


0.05 


0.05 


0.1 


retention 


rCDM 


0.2 


0.15 


0.03 


0.1 


ejection 



In Fig. 3 we show the present-day distribution of galax- 
ies in our ACDM model for this choice of parameters. Galax- 
ies are overplotted on the dark matter distribution using 
symbols whose size increases with luminosity and whose 
colour represents B — V index. The geometry of the slice is 
the same as in Fig. 1. Note the predominance of red galaxies 
in clusters and that of blue galaxies in the field. Notice also 
that the local voids are indeed nearly empty of galaxies. 



3. 5. 1 Tully-Fisher relations 

In K99 the Tully-Fisher relation was plotted for central 
galaxies with 'spiral morphology' assuming the disk rota- 
tion velocity Vdisk to be equal to V200 the circular velocity 
of the host halo. When we used the same prescription, we 
found a number of low luminosity galaxies, off by more than 
half a magnitude both from the observed relation and from 
the mean relation for the simulated galaxies. Similar out- 
liers can be seen in fig. 6 of K99 - at V200 ~ 300 km s~^ 
some of their 'Sb/Sc' galaxies are too faint by more than 
one magnitude. This problem arises because V200 ~Vdisk is 
a poor approximation for a subset of galaxies with unusual 
halo assembly histories. To get a better estimate for Vdisk, we 
proceed as follows. Throughout the evolution of our galaxy 
population we accumulate AMgtars V200 for each star-forming 
central disk where V200 is the circular velocity of its current 
halo. The circular velocity of the disk at any redshift is then 
defined by 

,.2 Z^, . ^-^Btars,z=z, ^2002=21 

Klisk = T7 , (6) 

JWstars,disk 

where the sum is over all redshift intervals since the last 
major merger and Mgtars.disk is the stellar mass of the disk. 
Satellite galaxies keep the disk rotation velocity they had 
when they last were central galaxies. 

In Fig. 4 we plot Tully-Fisher relations for simulated 
galaxies with 100 < Vdisk < 300km s^^ and 1.5 < Mb, bulge — 
Mb, total < 2.2 (appropriate for Sb/Sc galaxies according to 
Simien & dc Vaucoulcurs 1986). This is the same criterion 
used by K99 except that here we do not restrict ourselves to 
central galaxies. We compare our results with the fit to the 
observational data given by Giovanelli et al. (1997): 

M/-51og(/i) = -21.00±0.02-7.68±0.24(log(iy)-2.5),(7) 



where W = 2 x Vdisk. The dashed lines show zblo" of what 
Giovanelli et al. term the "intrinsic scatter" of their relor 
tion. This relation has been corrected for internal extinction 
so wc have not included any effects of dust in our model 
predictions. 

Even with our revised definition of Vdisk, there are still 
a few galaxies with low /-band luminosity which show up 
below the observed relation. We could remove the five low- 
luminosity galaxies at log(VI^) ~ 2.4 in ACDM by adopting a 
more stringent morphological selection criterion, for example 
1.6 < Mb , bulge MB.total < 2. 

We have checked that the galaxies lying below (more 
than 0.7 magnitude) the expected relation are redder than 
"typical" spirals: their mean B — V index is 0.9 and 1.0 
in ACDM and rCDM respectively, instead of the 0.7 found 
for "normal" spirals. These galaxies are older systems, for 
example, halo satellites which still satisfy the morphological 
selection criterion in the B-band, have kept the disk rotation 
velocity when they last were a central galaxy, but where star 
formation has stopped because of gas exhaustion, resulting 
in a relatively low stellar mass and 7-band luminosity. They 
might thus correspond to a subset of observed SO's. 

Giovanelli et al. (1997) quote an intrinsic scatter of 
tint = —0.28 X + 0.26 for their Tully-Fisher relation, with 
X = log W — 2.6. This is consistent with our ACDM model, 
but is exceeded by the tCDM model at the bright end. We 
find the dependence of the scatter on our rather uncertain 
morphological selection criterion to be quite strong. A larger 
scatter is expected in high-density cosmologies because of 
their late structure formation, as pointed out by Buchalter 
et al. (2001). 

3.5.2 Luminosity functions 

In Fig. 5 we compare the &j-band luminosity functions (LF) 
of our simulations with those derived from the 2dF and 
SDSS redshift surveys. The vertical lines correspond to our 
luminosity and morphology resolution limits. As already 
stressed, a major difficulty in reproducing observed lumi- 
nosity functions with the kind of modelling employed here 
is in establishing a sufficiently strong cut-off at high lumi- 
nosities by avoiding the formation of overly bright central 
galaxies in groups and clusters. Despite our best efforts and 
our ad hoc cooling switch at V200 = 350 km s^^, this has not 
been completely successfully achieved in the models, where 
we still have a bright galaxy at Mi^ ~ —24.5, possibly due 
to the overmerging issue discussed in Section 3.4. A second 
major difficulty is in suppressing the luminosity of galaxies 
in low mass haloes sufficiently to reproduce the fiat faint- 
cud slope of observed LFs. We have not stressed this issue 
here because our resolution does not allow us to address it 
adequately. 

The ACDM LF shows a deficit of L* galaxies, with re- 
spect to the observations, whereas the tCDM LF produces 
somewhat too many such galaxies. As pointed out by KaufT- 
mann et al. (1993), this difference is a consequence of the 
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Figure 4. TuUy— Fisher relation for Sb/Sc galaxies in the simulations. The solid line is the relation which Giovanelli et al. 
(1997) fit to their observational data. The dashed lines give their estimate for the "intrinsic scatter" of the data. 



differing number of V200 ~ 200 km s^^ haloes in the two cos- 
mologies, given our requirement that both should reproduce 
the observed Tully-Fisher relation. The global shape of the 
tCDM LF agrees somewhat better with the observations. 
The shapes of both our luminosity functions are closer to 
the observations than those given in K99. This is a result of 
our tuning our cooling and merging prescriptions as noted 
in previous sections (see also SP99 and COO) . 

The results of SOO suggest that the luminosity func- 
tion shape would be improved if we used simulations of 
sufficiently high resolution to follow the merging of galax- 
ies directly, rather than having to use a phenomenological 
model. The "knee" in our simulated &j-band LPs occurs at 

~ — 22 in both our models. 93% and 96% of the galax- 
ies brighter than these limits are central galaxies of haloes 
in ACDM and rCDM respectively. Of these central galax- 
ies, 73% and 16% belong to haloes with V200 > 350 km s~^, 
showing that our LF results in this range axe being affected 
strongly by the cooling cut-off and by the merging which 
can brighten galaxies above it. 



3.5.3 Masses for Milky Way look-alikes 

We select Milky Way look-alikes as the central galaxies of 
haloes with 200 <Vdisk < 240 km s~^ which also satisfy the 
"Sb/Sc" criterion we used when plotting the Tully-Fisher 
relation. Table 3 presents their mean stellar and cold gas 
masses together with their _B-band and /-band absolute 
magnitudes, their colours, their disk SFR's and their abun- 
dance in the simulations. 

Our two models show similar cold gas and stellar 
masses, quite comparable with those estimated for the Milky 
Way, although their present-day star formation rates differ 
by a factor of two. We find an average SFR of 1.8 M© /yr for 
disks of Milky Way-like galaxies in ACDM and 3.4 Mq /yr 
for such disks in rCDM. This reflects the later formation of 



structure in rCDM which requires more star formation at 
late times to get enough stars to reproduce the Tully-Fisher 
relation. Both values are in reasonable agreement with the 
2 to 3 Mq/ji estimated for the mean star formation rate 
in the Milky Way's disk over the last few Gyr (Rocha-Pinto 
et al. 2000ba). 

The difference in the number of galaxies with 200 <Vdisk 
< 240 km s"^ in ACDM and rCDM is substantial but is 
fully accounted for by the difference in the halo mass func- 
tions. 



3.5.4 Morphologies 

As in Kauffmann et al. (1993), we assign a morphologi- 
cal type based on the _B-band bulge-to-disk ratio: galaxies 
with //B.buige/ie.disk > 1.52 are ellipticals, SO's have 0.68 < 
ie.buige/ie.disk < 1.52, Spirals have iB.buige/i^B.disk < 0.68 
and irregulars have no bulge. 

We adjust our major merger threshold parameter so 
that the overall morphology fractions in our simulations 
match those derived by Baugh et al. (1996) from the observa- 
tions of Loveday (1996): 13%, 20% and 67% for Es, SOs, and 
Sps + Irrs. We estimate these ratios in the simulations using 
all galaxies brighter than our morphology resolution limit. 
Taking /bulge = 0.1 we obtain 12%, 19% and 69% as the rel- 
ative fraction of the three morphological classes in ACDM, 
and 16%, 14% and 70% in rCDM. For rCDM our SO frac- 
tion is somewhat smaller than the observed 20%. This may 
be another hint for overmerging: central galaxies generally 
spend little time within the SO bulge-to-disk ratio limits, 
either growing bigger disks or evolving into Es via merg- 
ers. Both the cooling and the merging processes are more 
rapid in rCDM, leading to a reduced rmmber of SOs. Given 
the different merging prescriptions used by different authors, 
one must be careful when comparing the /bulge values they 
derive. For example, SP99 and K99 both fit the observed 
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Figure 5. fej-band luminosity functions of the simulations. Dashed and dotted lines show the SDSS and 2dF GRS luminosity 
functions respectively. The vertical lines show the resolution limits in this band. The SDSS data have been converted to the 
2dF band as recommended by Blanton et al. (2001). 



Table 3. Masses (in units of Mq), disk star formation rates (in units of Mq/yi), colour and magnitudes of Milky- Way 
type galaxies (200 <VdiBk < 240 km s"!) 



Model 


Star mass 


Gas mass 


disk SFR 


Mb - 5 log h 


Mi - 5 log h 


B-V 


Number 


ACDM 


9.32 X 10^° 


1.05 X 10^° 


1.78 


-20.02 


-22.01 


0.74 


229 


tCDM 


1.58 X 10" 


1.15 X 10^° 


3.37 


-20.04 


-21.99 


0.71 


941 



morphological ratios with /bulge = 0.25, but with a differ- 
ent formula for the dynamical friction timoscale than is used 
here. As COO emphasize, the merger simulations of Walker 
et al. (1996) and Barnes (1998) suggest 0.1 < /bulge < 0.3. 

In Fig. 6, we plot mean morphology fractions as a func- 
tion of the maximum B absolute magnitude of the galax- 
ies considered. The fraction of bulgc-dominatod galaxies in- 
creases with luminosity. The most luminous galaxies, which 
are mainly the central galaxies of clusters, are all ellipti- 
cals. The variation of the morphological mix with luminos- 
ity is smooth, except for a steep rise in the elliptical frac- 
tion at Mb --^ —22. This is a consequence of our ad hoc 
suppression of cooling onto central galaxies in haloes with 
V200 > 350 km s^^. This suggests that the observed pre- 
dominance of bright ellipticals in clusters is related to the 
observed inability of cooling flows to generate new disks in 
such systems. 

3.5.5 Colour distributions 

Although we do not use colours when setting the free pa- 
rameters of our models, we briefly discuss here the colour 
distributions of our galaxies. The solid histogram in Fig. 7 
shows the distribution in B — V colour of all the galaxies 
brighter than our morphology limit. Two peaks are evident 
in these distributions which are separated when we plot in- 



dividual distributions for ellipticals and spirals, as defined 
in 3.5.4. There is good agreement with K99 and with the 
colours derived from the RC3 catalogue of de Vaucoulcurs 
et al. (1991) (shown as the dashed histogram). If wo restrict 
ourselves to galaxies residing in massive clusters (those with 
1^200 > 500kms~^), we obtain the dotted plots, showing 
that cluster galaxies are redder than their field counterparts 
(see also fig. 16 of SOO). 



4 NEARBY CLUSTERS 

This section begins the detailed comparison of our simula- 
tions with the individual structures seen in the local uni- 
verse. We start with a one-to-one comparison of the largest 
simulated clusters with their nearby counterparts. It is sur- 
prising that this comparison is possible at all because rich 
clusters arc systematically under-represented in the IRAS 
1.2 Jy survey from which we constructed our initial condi- 
tions. The next subsection shows that a convincing identifi- 
cation of simulated and real clusters is, nevertheless, possi- 
ble, and wc go on in subsequent subsections to explore this 
correspondance in terms of cluster dark matter properties, 
galaxy luminosity functions, and cluster mass-to-light ra- 
tios. 
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Figure 6. Morphology distributions for galaxies brighter than a given absolute Mg threshold. 
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Figure 7. Colour distributions for simulated galaxies brighter than the morphology limit. The left panels show all galaxies, the central 
panels ellipticals, and the right panels spirals. The top plots arc for ACDM, the bottom ones for rCDM. Solid histograms refer to the full 
simulated volume, dotted histograms to galaxies in clusters and the dashed histogram to the RC3. 
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4.1 Observed versus simulated clusters 

Among nearby rich clusters we select Coma, Virgo, Cen- 
taurus, Hydra and Perseus as particularly interesting for 
comparison with our simulations. All are relatively massive 
and lie within 8000 km s~^. We identify them with their 
simulated counterparts as follows. We take the observed 
coordinates and redshift from the nasa/ipac extragalactic 
database ^, and then look in the simulations for the most 
massive group less than 7 Mpc away from the expected 
redshift-space position. This search radius has been chosen 
arbitrarily to ensure that we recover all five selected clusters. 

The positional agreement between real and simulated 
clusters, given in Table 4, is quite good except for our ACDM 
Contaurus cluster, which is 6.7 Mpc away from the ex- 
pected position. The very good match for Coma is particu- 
larly surprising: the error in position is less than 2 Mpc 
in both models, despite the fact that the cluster is close to 
the edge of the high-resolution region, where the constraints 
arc noisy due to the relatively poor sampling of the 1.2 Jy 
survey. 

4.2 Cluster dark matter 

In Table 4 we list for each cluster the total dark matter 
mass (as given by our group- finder) , the virial mass (de- 
fined as the mass within ii20o), both in units 10^* h^^ Mq, 
the virial radius in Mpc, and the line-of-sight galaxy 
velocity dispersion in units of km s~^. We take the real data 
from Girardi et al. (1998). In the observed virial mass col- 
umn, the first figure is M2Q0 , inferred from the virial theorem 
applied to galaxies within 7i2oo but without including a sur- 
face term (the traditional virial estimate). The second figure 
is M2oo,corr, wMch iucludcs a term to correct for the surface 
"pressure" and should be more accurate (see Girardi et al. 
1998 for details). We compute lino-of-sight velocity disper- 
sions in our simulations using all galaxies above our resolu- 
tion threshold. The first figure is obtained by considering all 
galaxies in the cluster halo, the second only those galaxies 
within the virial radius. There is no significant difi^erence. 

It is striking that the ranking of clusters by mass in the 
simulations agrees with their ranking in the real universe. 
For individual clusters the simulated M200 is usually within 
a factor of 2 of the observed value and is often better. If we 
assume that M2oo,corr is the better observational estimator, 
then for ACDM Coma, Hydra and Perseus agree well, but 
Virgo is too massive by a factor of two, while for tCDM 
Virgo agrees well but Perseus and Coma are too massive by 
a factor of two whereas Hydra is undermassivc by a factor 
exceeding 2. This casts some doubt on the identification of 
Hydra in tCDM - however, it lies only 4 h^^ Mpc away 
from the expected location and the next major (and more 
massive) cluster is ~ 7 Mpc away. 

For the line of sight velocity dispersion the ranking of 
all five clusters agrees with observation in both cosmologies. 
The individual values also agree quite well in most cases. 
The larger deviations correspond to the discrepant masses 
just discussed. Recall that Centaurus has a well-known bi- 
modal velocity structure (Lucey et al. 1986), reflecting the 

t http://nedwww.ipac.caltech.edu/forms/byname.html 



Table 5. Apparent magnitudes and colours of brightest cluster 
galaxies. The first line is ACDM, the second one tCDM and the 
last one the observations. 



Name 


rriB 


mv 


B-V 


Coma 


11.57 

11.79 

12J>{) 


10.57 

in. 81 

11. Til 


1.00 

0.98 

().!)9 


Virgo 


9.64 
8.94 
9.42 


8.65 
7.94 
8.45 


0.99 
1.00 
0.97 


Centaurus 


10.47 
10.68 


9.48 
9.68 
10.12 


0.99 
1.00 


Hydra 


11.17 
12.42 


10.18 
11.45 
10.84 


0.99 
0.97 


Perseus 


10.61 
10.79 
11.79 


9.62 
9.80 
11.27 


0.99 
0.99 
0.52 



" AbcU ct al. (1989) 
Sandage & Hardy (1973) 

fact that the groups Cen30 and Ccn45 are currently merg- 
ing (Churazov et al. 1999). The quoted observed value is 
inferred if the two groups are considered together. 

4.3 Luminosities and M/L ratios 

For each cluster, we identified the most luminous galaxy (in 
B) and calculated its B — V colour. In Table 5, we compare 
the results with the data of Sandage (1972, hereafter S72). 
For each cluster the first line refers to the brightest clus- 
ter galaxy in our ACDM model, the second to the BCG in 
rCDM, and the last to the observed BCG. In (almost) all 
cases our simulated BCGs are brighter than observed. (The 
BCG in the low mass tCDM 'Hydra' is an exception.) The 
discrepancy is less than about a magnitude for all objects 
except the well-known peculiar galaxy NGC 1275 in Perseus, 
which is more than one and a half magnitude fainter than 
our prediction. The observed values quoted have been cor- 
rected for galactic extinction as follows. Most corrections 
computed by S72 used a simple, latitude dependent model, 
except for Perseus where S72 estimated the correction based 
on the B-V index of NGC 1275, and found Ay = 0.3 
and Ab = 0.4. Given the peculiar spectrum of the galaxy, 
these values are quite uncertain. The V-band magnitude of 
the BCG of Centaurus found in Abell ct al. (1989) has not 
been corrected for extinction. Recall that, owing to their low 
galactic latitude, extinction is a major issue for the B-band 
magnitude of Centaurus and Perseus. For these two clusters, 
we have therefore taken the estimation of galactic extinction 
given by Schlcgcl et al. (1998) from their 100 /im all-sky map 
of dust emission: they find Ay = 0.38 and Ab = 0.49 for 
Centaurus and Ay = 0.57, Ab = 0.74 for Perseus. 

The discrepancy in brightness between the simulated 
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Table 4. Dark matter in nearby clusters. The first line is ACDM, the second tCDM and the last the real data. 



Name 


SGL 


SGB 


cz 




M200 


.R200 




Coma 


90.2 


7.10 


7161 


6.56 


5.33 


1.32 


905,972 




90.0 


8.65 


7137 


14.8 


11.7 


1.72 


1202,1261 




oy .o 


o.oz 








1 64 


821 


Virgo 


107 


-11.4 


1413 


4.82 


4.01 


1.20 


790,820 




89.0 


13.1 


1221 


2.69 


2.55 


1.03 


716.729 








"1 ("V7<1 






1 0(; 


1 ) . ) z 


Centaurus 


159 


-4.5 


3873 


10.4 


4.85 


1.28 


982,956 




158 


-8.0 


3539 


3.59 


3.50 


1.15 


846,855 




156 


-11.4 


3298 








791 


Hydra 


142 


-36.0 


3417 


4.23 


2.59 


1.04 


662,709 




132 


-39.7 


3796 


0.81 


0.77 


0.69 


490,497 




139 


-37.5 


3418 




2.8,1.9 


1.22 


610 


Perseus 


348 


-11.8 


4914 


15.2 


8.74 


1.56 


966,1036 




353 


-16.5 


5346 


26.0 


15.9 


1.9 


1293,1347 




348 


-14.1 


5486 




15.6,9.1 


2.1 


1026 



BCGs and the data may be due in part to the luminosi- 
ties of the observed BCGs being underestimated (c.f. Uson 
et al. 1991; Gonzalez et al. 2000), and in part to the over- 
merging problem discussed in 3.4. Note that Coma has two 
very bright galaxies in a binary in its core. The luminosities 
of our simulated Coma BCGs are similar to the sum of those 
of NGC 4889 and NGC 4874. 

The colour indices of our simulated BCGs arc all very 
red, reflecting their old stellar populations. (Note that we 
use solar metallicity population synthesis models when esti- 
mating the observed properties of our stellar populations.) 
They agree well with the observed BCG colours in all cases 
except that of NGC 1275. The anomalous A-typo spectrum 
of this galaxy is thought to result from star formation asso- 
ciated either with the strong cooling flow in Perseus or with 
an apparent ongoing merger with a gas-rich galaxy. 

In Fig. 8, we give my for the ten most luminous galaxies 
in our simulated Coma and Virgo clusters. They are com- 
pared to the data of Gudehus (1995). The BCG has been dis- 
cussed in the previous section, and is typically brighter in the 
models than in the simulations. The remaining nine galax- 
ies are fainter than observed by about one magnitude. It 
appears that choosing parameters to fit the TuUy-Fisher re- 
lation has resulted in too few stars forming in cluster galax- 
ies. Possible reasons might be a systematic underestimate 
of Vdisk (and thus the assignment of low L) for the central 
galaxies of our haloes, or an underestimate of the effects of 
dust (and thus of the stellar masses) in these same 'TuUy- 
Fisher' galaxies. The first explanation would affect the lumi- 
nosities of all galaxies and would so shift the overall LFs of 
the models, the second would affect primarily cluster galax- 
ies and would have little effect on the global LFs. The fact 
that Diaferio et al. (2001) get a good fit (actually slightly 
overluminous) to the bright end of the LF's of the CNOCl 
clusters, while SOO slightly underpredict the total luminosity 
of Coma, suggests that the problem is not fundamental but 
lies in the details of the phenomenological modelling. The 
discrepancy in our own models is also evident in the cluster 
mass-to-light ratios as we now show. 



Table 6. Mass-to-light ratios T of clusters. 



ACDAl 



rCDAl 



Data 



Name 


Ti 


T2 


Ti 


T2 


T2 


Coma 


514 


614 


733 


762 


121-225 


Virgo 


485 


558 


770 


824 


249-459 


Centaurus 


470 


574 


875 


892 




Hydra 


484 


657 


629 


651 


162-327 


Perseus 


376 


377 


578 


703 


280-458 



We compare simulated mass to light ratios in the bj- 
band with those derived by Girardi et al. (2000). Wc ap- 
proximate the bj magnitude by taking the "galaxy colour 
relation": bj = B - 0.35 {B - V) (Blanton et al. 2001). We 
estimate two ratios for our simulated clusters, one using the 
total mass and light in our simulated haloes, the other using 
values within their virial radii: 

Ti = (8) 



T2 



Lbj ,tot ' 

Mvir 



(9) 



We give all ratios in solar units for h = 1. 

If (pbj is our luminosity function for the &j-band, we 
include an estimated contribution to the total luminosity 
from galaxies fainter than our luminosity resolution limit as 
follows: 



Lbj ,tot — Lb J ,faint + Lb^ .bright 

ibj, bright = ^ ^ Lbj ,i 

Lfj . > resolution 



^6,-, faint 



Nbri 



ght 



<t>b,{L)dL 



Jo 



(10) 
(11) 

L4>bAL)dL (12) 
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Figure 8. Apparent magnitude in the V-band of tlie briglitest ten galaxies in the simulated Coma and Virgo clusters 
compared to observations. 



Table 6 lists the results: the two numbers cited for the data 
are the minimum and maximum values found by Girardi 
et al. (2000) when varying their sample of galaxies, their 
maximum clustercentric radius, and their method to esti- 
mate luminosities. For Centaurus, Girardi et al. (2000) did 
not provide any data since the cluster is multipeaked. The 
observed ratios of Girardi et al. (2000) in the 6j-band scatter 
around a value of 250. 

For our ACDM model, we find values scattering around 
Ti ~ 500 and T2 ~ 600, and around Ti ~ 600 and T2 ~ 
700 for the tCDM model, although with a noticeably larger 
scatter. 

The higher value within the virial radius is a conse- 
quence of the radial variation of the morphology and colour 
of the galaxies. The typical discrepancy with observation ap- 
pears to be factor of 2 and 2-3 for the ACDM and rCDM 
models respectively. This is somewhat larger than expected 
from the results presented above since our simulated cluster 
masses are not systematically in error, the brighter galax- 
ies in our clusters are typically too faint by a magnitude or 
less, and our cluster luminosity functions (pb^ have faint end 
slopes (which are important for the luminosity correction 
discussed above) which are steeper than those of most ob- 
servational determinations (Folkes et al. 1999; Blanton et al. 
2001). 

Further work, including careful analysis of the obser- 
vational mass determination techniques is needed to under- 
stand how much of this discrepancy comes from incorrect 
galaxy modelling and how much from systcmatics intro- 
duced by the mass and luminosity measurement techniques. 
For example, the higher M/L ratios gencrically obtained in 
the tCDM model compared to ACDM are not reflected in 
the Johnson _B-band, where the figures are similar. At these 
wavelengths, our ACDM simulation gives total M/L ratios 
of 536, 501, 512, 490 and 402 for Coma, Virgo, Hydra, Cen- 
taurus and Perseus respectively. Using modelling techniques 
very similar to ours, K99 found the S-band M/L ratios of 



clusters in ACDM to be biased with respect to that of the 
Universe as a whole by an amount which is quite consistent 
with observation. In an independent analysis Benson et al. 
(2000b) obtained Tb ~ 400 for massive haloes in ACDM. 
For comparison, Kent & Gunn (1982) measured Tb = 360 
for the Coma cluster, perhaps the best observed of local 
clusters. 

It is also possible to use the galaxy populations in 
our simulations to study the radial distribution of galaxies 
within clusters. Detailed analyses of this kind were carried 
out by Diafcrio et al. (1999 2001) comparing their GIF sim- 
ulations to CfA and CNOC clusters. At the resolution of our 
(and their) simulations, it is only possible to get useful re- 
sults by stacking many clusters. We have done this with our 
data and find results which are very similar to that of the 
earlier work. In the mean the radial distribution of galaxies 
within clusters is similar to but slightly less concentrated 
than that of the dark matter. Blue and star-forming galax- 
ies axe much less concentrated towards the cluster centres 
than redder galaxies. The proportion of elliptical galaxies 
rises (and that of spirals drops) near the cluster centre. We 
do not repeat the plots of the earlier paper here since we 
have little to add. Further progress on these issues requires 
simulations of significantly higher resolution, such as those 
of SOO. 



5 COMPARISON CATALOGUES 

In this section, we recall briefly the features of the two local 
galaxy catalogues with which we compare our simulations. 
We use the far-infrared selected PSCz catalogue since it in- 
cludes the IRAS 1.2 Jy survey which was used as the density 
constraint on our initial conditions. In addition it has near 
full-sky coverage and the highest available IR source den- 
sity in the region we have modelled. A complication is that 
comparison with this survey requires us to model the FIR 
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luminosity of simulated galaxies. At optical wavelengths, wc 
have chosen the recently completed UZC as our reference 
catalogue. Its sky coverage is substantially smaller than that 
of the PSCz but it has a higher source density and we can 
compare the observed B-band luminosities directly to our 
simulations. 



5.1 PSCz catalogue 

The PSCz catalogue has been described by Saunders et al. 
(2000). It contains 15411 IRAS galaxies and covers some 
84% of the sky, excluding regions of low galactic latitude and 
cirrus, and regions unobserved by the IRAS satellite. 1.2% 
of IRAS galaxies within the PSCz region and with bj < 
19.5 have unknown redshift, and we discard them from the 
catalogue. We then consider the PSCz catalogue is limited 
both at bj < 19.5 and at feoi^m > 0.6 Jy. Also, we only 
consider galaxies with cz < 8000 km s~^. 

To compare the results of our simulations with the 
PSCz, we need to estimate the IR luminosity of the galaxies 
at 60 /im. Our model as described so far does not provide us 
with thisg information. We assume that the 60 /im luminos- 
ity of a star-forming galaxy has two components, one com- 
ing from star-forming regions and proportional to its star 
formation rate, the other representing the re-emission of ob- 
scured light from older stars and proportional to its 7-band 
luminosity. For very low star formation rates we assume a 
negligible FIR luminosity. This threshold in star formation 
rate is set empirically. Formally, we write: 

M* 

Lfih = asFR M* + ai Li if M* > 0.01 (13) 

tHubble 

Lfir = otherwise (14) 

This is known to be at best a rough approximation for IRAS 
galaxies (Hclou 1986), but it gives sufBcicntly accurate re- 
sults for the purposes of this paper. A more detailed analy- 
sis would require consistent inclusion of the IR band in our 
SA model (see Devriendt & Guiderdoni 2000). We deter- 
mine the proportionality coefficients asfr and ai for each 
cosmology by matching the IR luminosity function of simu- 
lated galaxies brighter than the resolution limit to the IRAS 
luminosity function given in Saunders et al. (1990). Because 
of resolution effects we are able to fit the observed IRAS 
luminosity function only above a limiting 60 /itm luminosity 
of Leo = 2.7 IO^Lq and Leo = 1-50 lO^L© in ACDM and 
tCDM respectively, which we define as our FIR resolution 
limits. Here, Leo is expressed in units of the bolometric lumi- 
nosity of the sun, and is consistent with the assumed galaxy 
FIR spectrum of Saunders et al. (1990). To avoid confu- 
sion we will in the following refer to our previous luminosity 
resolution limits as the "optical resolution limits" to distin- 
guish them from the "FIR resolution limits" . Interestingly, 
the fitting values of asFR that we obtain: aspR = 1.68 lO'' 
and asFR = 6.4 10* in ACDM and tCDM respectively arc 
close to the one estimated by Rowan-Robinson (2000) for 
bright IR galaxies. As a check, we compare the number 
of simulated and observed galaxies with /eo/jm > 0.6 Jy, 
cz < 8000 km s~^, bj < 19.5 and brighter than both the op- 
tical and FIR resolution limits. The figures are reported in 
Table 7. Note that they are different for our two cosmologies 
because of the differing simulation resolution limits. 



5.2 UZC catalogue 

According to Falco ct al. (1999), the UZC is ~ 98% com- 
plete for mzwicky ^ 15.5 in the ranges 20*" ^ aigso ^ 4*^ 
and 8'' ^ 01950 ^ 17*', for declinations -2° .5 ^ Jigso ^ 50°. 
This covers roughly one third of the sky. We again consider 

only galaxies with cz < 8000km s~^, and with absolute mag- 
nitudes above the optical resolution limit for our simulated 
galaxies. 

To compare with the UZC, we need to cut the simu- 
lated galaxy catalogues in the appropriate S-band magni- 
tude. However, the assumption of Mzwicky ^ Mb is known 
to be rather approximate (Huchra 1976; Bothun & Cornell 
1990), with a 1 cr error of 0.3 mag. We adjust our threshold 
in simulated ms to obtain good agreement between the pre- 
dicted and observed number of galaxies in the UZC region. 
This process is independently applied to both our models. 
It leads us to adopt Mzwicky — Mb = —0.4 and 4-0.9 for 
the ACDM and rCDM models respectively. These shifts are 
larger than can plausibly be attributed to the difference in 
photometric systems, but allow us to compensate in part for 
the systematic offsets in luminosity function that are visible 
in Fig. 5. Again, the resulting number of galaxies in our real 
and mock UZC catalogues are reported in Table 7. 



6 SMOOTHED DENSITY FIELDS 

The constrained density field used to set up the initial con- 
ditions for our two simulations was derived from the IRAS 
1.2 Jy catalogue, smoothed with a Gaussian of dispersion 5 
Mpc. As a result we expect the galaxy density fields in 
our simulation to resemble each other and to resemble the 
real galaxy field when smoothed on similar scales. In this 
section we will check how well this works out. In the fol- 
lowing, we will denote by "optical galaxies" the sample of 
all simulated galaxies brighter than our S-band resolution 
limit. Similarly, we will call "FIR galaxies" all simulated 
galaxies brighter than both the optical and the FIR reso- 
lution limits. We will apply apparent magnitude limits to 
these samples only in section 6.4 where we compare with 
the observed PSCz sample. 

6.1 Smoothed maps 

As a qualitative example. Figs. 9 and 10 show the over- 
densities of optical galaxies in the supergalactic plane, us- 
ing Gaussian kernels with 5 and 10 Mpc smoothing 
length respectively (the symbol R in Equation 15 below). 
As expected, there is already good agreement between the 
two cosmologies on the 5 h^^ Mpc smoothing scale, and 
the match is even better at 10 Mpc. The three ma- 
jor overdensity peaks are Coma at the top of the plot, the 
Great Attractor at the center-left, and the Pisces-Perseus 
complex at the center-right. The Virgo cluster/Local Su- 
pcrclustor complex is also evident just above the centre of 
the plots. These maps can be compared with galaxy den- 
sity maps reconstructed directly from various IRAS surveys 
(e.g. Brancliini et al. 1999). The resemblance is quite close 
and is encouraging because of the many steps between using 
the observed density fields to constrain the initial dark mat- 
ter distributions in our models and constructing simulated 
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Table 7. Number of galaxies in mock and real PSCz and UZC surveys. 



Model 


mock PSCz catalogue 


PSCz data 


mock UZC catalogue 


UZC data 


ACDM 


6806 


6735 


8061 


8031 


rCDAl 


5 110 


5 112 


711)8 


7121 




galaxy density fields from the galaxies which are formed dur- 
ing their evolution. 

To obtain more quantitative results, we now compare 
various smoothed density fields point-by-point. To do this we 
assign the dark matter particles and the galaxies to a regular 
grid with a spacing of 10 Mpc using a CIC scheme. 
We keep only those grid points within cz < 8000 km s~^, 
resulting in some ~ 2100 nodes. We further smooth this 
density field using a Gaussian kernel of the form 



(15) 



We take R = 10h~^ Mpc, corresponding to twice the 
smoothing used when generating the initial conditions. Af- 
ter this smoothing, we correct for artificial effects introduced 
by the spherical boundary of our galaxy distribution. This 
is done by multiplying the smoothed field by the ratio of 
fields obtained by applying similar smoothing procedures to 
two Poisson distributed catalogues, one uniform within a 
box of side 220 Mpc centred on the simulation sphere, 
the second uniform within the region where our data lie 
{cz < 8000 km s~^, together with the PSCz mask when ap- 
propriate) . We also evaluate the volume of the intersection 
between the smoothing kernel and our selection region using 
this simple Monte-Carlo scheme and we discard nodes where 
this volume is less than half of the volume of the smoothing 
kernel. 

Note that we have chosen the spacing of our grid so 
that the smoothed density fields are not heavily oversam- 
pled, while retaining a reasonable number of data points. For 



consistency, we define the mean density of each smoothed 
field as the straight average of its values at our final set of 
nodes. We can then define the smoothed overdensity at each 

node through 



ps.lO 



Ps.lO 



(16) 



We now discuss scatter plots which show point-by-point 
comparisons of various of these density fields. We begin by 
focussing on difTerent components within a given simulation. 
We then compare each component between our two simula- 
tions. Finally we compare our simulated FIR galaxy density 
fields with the real PSCz data. 



6.2 Optical galaxies vs. FIR galaxies vs. mass 

In Fig. 11, we compare the overdensity distributions of op- 
tical galaxies, of FIR galaxies and of mass in each of our 
simulations. 

The overdensities of optical galaxies and of dark mat- 
ter are remarkably tightly correlated in both simulations. 
In each case the galaxies are almost unbiased with respect 
to the underlying DM except in the highest density regions 
(5s, 10 ~ 1 to 3), where the optical ACDM galaxies are sig- 
nificantly antibiased. 

Biases are more evident and the scatter is larger when 
the FIR overdensity fields arc compared with those for 
the dark matter. For ACDM the trends are similar to but 
stronger than those for the optical galaxies. At DM overden- 
sities of order (5s, lo ~ 3, galaxies are antibiased by a factor of 



© 0000 RAS, MNRAS 000, 1-24 



16 H. Mathis et al. 




almost 3. This is expected given our model for the FIR lu- 
minosity: cluster galeixies have low star formation rates and 
so are assigned small or zero Lfir (consistent with observa- 
tions of most E/SO galaxies) and so arc almost all excluded 
from our simulated FIR catalogues. This substantially re- 
duces the number of galaxies counted in the densest regions. 
On the other hand, in mean density and underdense regions 
star-forming spirals predominate, and these are included in 
both the optical and the FIR samples. As a result the scat- 
ter plots for the two populations are quite similar in this 
density range. 

A similar but weaker pattern is observed for the rCDM 
FIR galaxies and for the same reasons. The weakening in 
this case reflects the fact that the present day star formation 
rate in the tCDM cosmology is five times that in the ACDM 
cosmology, and a significant number of galaxies in dense 
regions have fallen in recently enough for their ongoing star 
formation to take them above our 60 /um resolution limit. 

The larger scatter in the FIR plots may well be a real 
effect, reflecting, perhaps, the stochasticity introduced by 
the fact that FIR luminosity is linked to star formation rate 
rather than stellar mass. The resulting underweighting of all 
group and cluster environments introduces scatter as well as 
bias into the relation between smoothed mass and smoothed 
light. We checked that the lower scatter in the optical case is 
not a result of the larger number of galaxies in our simulated 
optical samples. Randomly picking one third of the optical 
galaxies and using this subset to construct the overdensity 
field does not lead to a noticeable increase in the scatter in 
the optical plot. 

As a consistency check, we plot in the two bottom 
panels of Fig. 11 the overdensities of FIR galaxies against 
those of optical galaxies. The effect of our FIR modelling 
is again very clear: in both models, but particularly in the 
ACDM model, significantly overdense regions in terms of op- 
tical galaxies become moderately overdense in terms of FIR 
galaxies. 



6.3 ACDM vs. rCDM 

In this section, we compare the distributions of dark matter, 
of optical and of FIR galaxies between our two cosmologies. 
Recall that although the initial conditions of the two sim- 
ulations are nearly identical when smoothed on scales of 5 

Mpc (apart from a difference in fluctuation amplitude) 
they differ on smaller scales. As a result, the build-up of 
galaxies is almost uncorrelated between them. 

The upper left plot of Fig. 12 compares dark matter 
overdensities in the two models. The shape is as expected 
from the initial set-up of the simulations. Initial amplitudes 
were chosen so that the abundance of objects of rich clus- 
ter mass would agree in the two cases. This requires the 
amplitude of linear overdensity fluctuations to be larger in 
ACDM than in rCDM. The difference in amplitude carries 
over to the flnal time and shows up as a clear "antibias" of 
the rCDM density relative to the ACDM field. On the other 
hand, the correlation between the two overdensity fields is 
strong, showing that the differences in the initial conditions 
(and so in the evolved mass fields) on small scales have lit- 
tle influence on larger scale fluctuations for the 10 Mpc 
smoothing used to make these plots. 

The upper right plot shows that the galaxy formation 
recipes we have implemented, which are tuned to reproduce 
the observed TuUy-Fisher relation and, to a lesser extent the 
observed luminosity function, are only partially successful in 
compensating for this difference in dark matter clustering. 
The distribution of optical galaxies in rCDM is still signif- 
icantly antibiased relative to that in ACDM, although the 
effect is weaker than for the dark matter. 

Somewhat surprisingly, this difference in clustering am- 
plitude is almost absent in the simulated FIR galaxy distri- 
butions. The lower left plot in Fig. 12 shows that the over- 
density fields of the two simulations are very similar for this 
population. The higher star-formation rates of non-cluster 
galaxies in the rCDM model clearly boost the relative am- 
plitude of the density fiuctuations in FIR population sufli- 
ciently to compensate for the weaker clustering of galaxies 
by optical luminosity (or stellar mass) . 



© 0000 RAS, MNRAS 000, 1-24 



Simulating the Formation of the Local Galaxy Population 




ACDM dark matter xCDM dark matter 




ACDM dark matter xCDM dark matter 




ACDM optical Galaxies xCDM optical Galaxies 

Figure 11. Scatter plot comparisons of the overdciisity fields for optical galaxies, for FIR galaxies and for mass in 
our two simulations. The top row plots optical galaxies against mass, the middle row FIR galaxies against mass, 
and the bottom row FIR galaxies against optical galaxies. ACDM is on the left and tCDM is on the right. The 
diagonal line is y = a; in each plot. 
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In the lower right plot, we compare the redshift-space 
distributions for this same FIR galaxy population. As ex- 
pected, the scatter is greater in redshift space and the over- 
all amplitude of the density fluctuations is increased. In ad- 
dition, an extended plume appears at high densities. These 
points correspond to the environments of massive clusters, in 
particular the Great Attractor and Pisces-Perseus regions. It 
is interesting that peculiar velocities appear to enhance the 
density contrast of these structures more strongly in ACDM 
than in rCDM. 

From this series of four plots we conclude that the point- 
by-point agreement between our two simulations is quite 
good for all three components, at least for the rather large 
smoothing employed here. As expected, the mass is more 
weakly clustered in rCDM, and this effect carries over to 
the optical galaxy population. The smoothed FIR galaxy 
populations in the two simulations are very similar. 



6.4 Simulated PSCz galaxies vs. observations 

We now compare the distribution of the simulated FIR 
galaxies directly to the observed distribution in the PSCz 
catalogue. This requires a slightly different FIR sample than 
that used in the last section. In addition to the redshift and 
simulation resolution limits already enforced, we need to ap- 
ply the sky mask of the observed catalogue and its flux limits 
at 0.6 Jy and hj < 19.5. The effects of these additional selec- 
tion criteria are consistently taken into account when con- 
structing overdensity fields both for the observations and for 
the simulated data. 

The left and right plots of Fig. 13 compare simulated 
and observed overdensity fields for the ACDM and rCDM 
cosmologies respectively. These fields are, of necessity, con- 
structed in redshift-space. The agreement is reasonably 
good, with greater scatter in the ACDM case. In neither 
cosmology is there any obvious bias between simulation and 
observation. This is a nice confirmation that our SA and FIR 
luminosity modelling schemes produce an FIR galaxy pop- 
ulation with a large-scale distribution which quite closely 
resembles that used to define our initial conditions. 



7 CORRELATION FUNCTIONS 

In this section we use correlation techniques to extend both 
our clustering measurements and our comparison with ob- 
servations to smaller scales. We begin by using auto- and 
cross-correlations of the various simulated populations to 
quantify small-scale bias. We then cross-correlate observed 
and mock catalogues to further explore how well we have 
reproduced our own cosmic neighborhood. 



7.1 Autocorrelations 

In Fig. 14 we plot real-space autocorrelation functions for 
the dark matter and for galaxies brighter than our B-band 
resolution limit ( "optical galaxies" ) for each of our simula- 
tions. For reference, we also plot the dark matter autocor- 
relations for the GIF simulations analysed by K99 and the 
autocorrelation function of real optical galaxies as inferred 
from inversion of the angular correlation data for the APM 
survey (Baugh 1996). We compute our correlation functions 



up to a scale of 10 Mpc. Note that the analysis in Jenk- 
ins et al. (1998) shows that the GIF results are close to the 
ensemble-averaged expectations for the two cosmologies. 

A comparison of the dark matter correlations between 
our simulations and the GIF simulations shows similar be- 
haviour in the two cosmologies. There is reasonably good 
agreement on small scales, but our simulations have more 
power than the GIF models on scales of a few Mpc. This is 
clearly a reflection of the particular way in which we have 
constrained the large-scale density field in our initial condi- 
tions using the observed distribution of 1.2 Jy galEixies. 

On small scales the optical galaxy distribution is quite 
strongly antibiascd relative to the mass in both our cosmolo- 
gies, while on scales of a few Mpc the difference is much 
smaller. Indeed, in tCDM the correlations of dark matter 
and optical galaxies are almost equal beyond 1.5 Mpc. 
For neither case are the autocorrelations of the optical galax- 
ies close either to a power law or to the APM data. It is un- 
clear whether this is a problem, since correlations calculated 
for the relatively local region we are modelling should not 
necessarily reproduce those found for much larger "represen- 
tative" regions. We will see below that our ACDM model 
docs seem to reproduce the correlation statistics of local 
galaxies quite well, while our rCDM model does not. On 
scales of 5 to 10 Mpc the correlation amplitude for op- 
tical galaxies is quite close to the APM values in both our 
simulations. 



7.2 Cross-correlations 

To make a more quantitative comparison between our simu- 
lations and the observed PSGz and UZC catalogues we have 
made mock PSCz and UZC catalogues which reproduce in 
detail the sky coverage and the apparent luminosity limits 
of the observational data. We also limit both the mock cat- 
alogues and the real catalogues to galaxies with absolute 
luminosities brighter than the relevant resolution limits of 
our simulated catalogues (see above). We can then compare 
real and simulated distributions in detail using auto- and 
cross-correlations. We have applied the same sky masks and 
depth selection functions to our simulated mass distributions 
to produce mass catalogues which can be compared directly 
with the observed galaxy distributions using the same tech- 
niques. 

To compute auto- and cross-correlations we use the es- 
timator suggested by Hamilton (1993): 

{DiD2){RiR2} ^ 

where (DD), (RR) and (DR) refer to the number of data- 
data, random-random, and data-random pairs, and the sub- 
scripts refer to the two catalogues. We compute these cor- 
relations from 0.7 Mpc up to 15 Mpc. The random 
catalogues used here are generated using the same angular 
mask as the corresponding galaxy (or mass) catalogue and 
a selection function in depth derived from the relevant lumi- 
nosity function (that of Saunders et al. 1990 for the PSCz, 
mock PSCz and associated mass catalogues; that of Marzke 
et al. 1994 for the UZC and its associated mass catalogue; 
those of 3.5.2 for the mock UZC catalogues). We ensure 
that each random catalogue contains at least ten times as 
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Figure 12. Scatter plots of the overdensities in rCDM against those in ACDM. The upper left plot shows the 
dark matter, upper right and lower left the optical and FIR galaxies in real-space, lower right the FIR galaxies 
in redshift-space. 



many points as its corresponding "data" catalogue. This en- 
sures that uncertainties in pair counts are dominated by the 
number of available galaxies (or mass particles). 



7.2.1 Cross-correlations with FIR galaxies 

In Fig. 15, we give redshift-space auto- and cross-correlation 
functions for our observed PSGz, mock PSCz and PSCz 
mass catalogues. 

A surprising result in this figure is that the autocorre- 
lation of the mock PSCz catalogue agrees almost perfectly 
with that of the real data for ACDM and is quite close to it 
for separations above 1.5 Mpc in the rCDM case also. 
In both cases the model autocorrelation is close to a power- 
law of index —1.28, the value found by Fisher et al. (1994) 
for the 1.2 Jy sample as a whole. Over the range we plot, 
the bias of the simulated PSCz galaxies relative to the mass 



is almost constant; the ACDM galaxies are substantially an- 
tibiased while the rCDM galaxies are almost unbiased. 

As expected, cross-correlations of the observed PSCz 
galaxies with either our mock PSCz galaxies or the simu- 
lated mass distribution show a different behaviour. On scales 
larger than about 5 Mpc the cross-correlation between 
real and mock galaxies is as strong as (and effectively equal 
to) the autocorrelation of either population. This is a strik- 
ing confirmation of the effectiveness of our constraint and 
galaxy formation procedures. On scales below 5 Mpc 
this cross-correlation flattens out to a constant value, reflect- 
ing the fact that the small-scale structure in our simulations 
is unrelated to that in the real Universe. 

The cross-correlation between the observed PSCz 
galaxies and the simulated mass distribution has a similar 
shape, and indeed, for tCDM (where the mock PSCz galax- 
ies are almost unbiased) it is almost identical to the cross- 
correlation with the mock galaxy catalogue. For ACDM the 
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Figure 13. Scatter plots of the redshift-space overdensities of simulated PSCz galaxies against those of the real 
PSCz data. 




Figure 14. Autocorrelation functions for dark matter and for galaxies brighter than our _B-band resolution limit. Solid, 
dotted, and dash— dotted lines are for the CR dark matter, the CR ojrtieal galaxies, and the GIF dark matter respectively. 
The error bars are the results obtained from angular correlations in the APM catalogue. ACDM is on the left and tCDM 
on the right. 



cross-correlation with mass is stronger than that with the 
mock galaxies. On large scale this cross-correlation is ap- 
proximately the geometric mean of the autocorrelations of 
the two populations. This would be expected for a pure lin- 
ear bias model but may be surprising in view of the scatter 
and tjhe nonlinearity evident in Figs. 11 and 13. 



7.2.2 Cross-correlations with optical galaxies 

Fig. 16 shows auto- and cross-correlation functions in the 
same format as Fig. 15 but for observed and mock galaxy 



catalogues and for mass catalogues with the sky mask and 
depth distribution of the UZC catalogue. 

It is interesting to compare the autocorrelation func- 
tions of this plot with the real-space autocorrelations in 
Fig. 14. As was also the case for the mock FIR galaxies, 
the redshift-space autocorrelation functions of the mock op- 
tical galaxies are parallel to those of the dark matter; for 
ACDM the galaxies are antibiased, although less strongly so 
than the FIR galaxies, while for rCDM the galaxies show 
a slight positive bias. In both cases the dark matter and 
galaxy curves are much more nearly parallel than are the 
corresponding real-space functions in Fig. 14. Note also the 
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Figure 15. Auto- and cross-correlation functions for tlie observed and simulated PSCz catalogues and for the mass. Note 
that the mass catalogues used here have the same sky and depth distribution as the observed catalogue. 



substantial differences between the mass autocorrelations 
plotted in Fig. 16 and those shown for the same simula- 
tions in Fig. 15. These reflect the different sky coverage and 
depth distribution of our PSCz and UZC catalogues and 
emphasise that neither should be considered a fair sample 
of the Universe as a whole. 

When wo compare our mock catalogue autocorrelations 
with those of the real UZC galaxies we find excellent agree- 
ment below 5 Mpc for ACDM with a slight underpredic- 
tion of the observed amplitude on larger scales. For tCDM 
the predicted autocorrelations arc low on all scales, with 
the difference in ^(s) ranging from 30% on large scales to 
a factor of 2 around 1 Mpc. Comparing with Fig. 14 
and Fig. 15, we see that while the autocorrelations of ob- 
served optical galaxies in our UZC region are significantly 
weaker than those of observed FIR galaxies in our PSCz re- 
gion, they nevertheless substantially exceed those measured 
for galaxies in the APM survey (even allowing for the differ- 
ence between real and redshift space) . This again emphasises 
that the regions of space we are analysing are quite small 
and should not be thought representative. 

The cross-correlations of the observed UZC galaxies 
with the mass and mock UZC galaxy distributions in our 
simulations are almost flat below ~ 5 Mpc, especially 
in the rCDM model. On larger scales they do approach the 
autocorrelation amplitude of the observed galaxies but the 
convergence is less compelling than was the case when we 
compared our mock FIR galaxies with the PSCz data. This 
difference may in part reflect the smaller volume covered by 
our UZC catalogues. It may also be due to the fact that 
we used an FIR galaxy catalogue when setting up our con- 
strained initial conditions. 

To summarise these sections on correlation functions, 
our initial condition generation and galaxy evolution proce- 
dures have allowed us to build physically consistent ACDM 
and rCDM models which reproduce well both the individual 
structures and the statistical clustering properties of star- 



forming galaxies in our local neighborhood. For ACDM the 
statistical properties of optically selected galaxies from the 
UZC catalogue arc also well reproduced, but this is not the 
case in our tCDM simulation where the optical galaxies are 
significantly more weakly clustered. 



8 MOCK CATALOGUES OF PECULIAR 
VELOCITIES 

As an example of how our simulations can be used to test 
and calibrate the methods of cosmic-flow analysis, we have 
generated Monte-Carlo mock samples mimicking real cat- 
alogues of peculiar velocities. In particular, we show here 
mock catalogues corresponding to the Mark III Catalogue 
of Peculiar Velocities (WiUick ct al. 1995 1996 1997). These 
were made with an improved version of the procedure de- 
scribed in K96, which provided the basis for earlier testing 
of methods. The major advance allowed by our new simu- 
lations comes in part from their ability to follow the non- 
linear evolution of the mass distribution on sub-Mpc scales, 
but more importantly from the fact that they simulate di- 
rectly the formation of the galaxies themselves. In K96, 
"galaxies" were inserted by applying a statistical biasing 
scheme to the smoothed mass density field, and their mag- 
nitudes were drawn at random from an assumed luminosity 
function. Here, galaxies form in a physically consistent way 
through the condensation of gas at the centres of dark haloes 
and their morphologies reflect their actual merging histories. 
This more realistic treatment should better capture possible 
correlations between galaxy properties and the underlying 
density and velocity fields. Our simulations, based on con- 
strained realizations, then allow a full reproduction of the 
correlation between systematic errors and the actual signal 
of density and velocity fields. 

For the purpose of this section we distinguished between 
"spiral" and "elliptical" galaxies by defining ellipticals as 
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Figure 16. Auto- and cross-correlation functions for observed and simulated UZC galaxies and for the simulated mass 
distribution. ACDM is on the left, rCDM on the right. 



galaxies with total stellar mass > 4.5 x 10 M©, and bulge- 
to-total V^-band luminosity ratio > 0.4. Small adjustments 
wore made to the semi-analytic absolute magnitudes in or- 
der to fine-tune their match with a Schechter luminosity 
function. We then assigned to each galaxy an "observed" 
linewidth based on an assumed TuUy Fisher relation and 
scatter as in K96 - in practice these liriewidths were al- 
ways close to those assigned by the algorithms described 
earlier which are based on the actual circular velocities of 
the haloes. To simulate the selection based on magnitude in 
each dataset of the Mark III catalogue, we assumed Galac- 
tic extinction as a function of Galactic latitude as in K96, 
and we slightly smeared the magnitude limits to take into 
account scatter in the relation between the magnitudes used 
for selection and those appearing in the TuUy-Fisher rela- 
tion. Fig. 17 demonstrates our success in matching the final 
distribution of redshift and apparent magnitude in our mock 
catalogues with that of the real data from the largest single 
dataset within the Mark III catalogue. 

Each mock dataset was then diluted at random (trying 
to mimic selection by other independent properties such as 
inclination or lack of a strong bar) to match the number 
of galaxies in the real catalogues. This random sampling, 
along with the random distance errors (introduced by the 
TF scatter) was repeated 10 times to generate 10 mock cat- 
alogues. The mock data were then grouped, and statistically 
corrected for Malmquist bias in exactly the same way as the 
real data (details are in the Mark III papers and are sum- 
marised in the POTENT analysis, Dekel et al. 1999). 

Mock and real maps of radial peculiar velocities in a 
slice about the Supergalactic plane arEe shown in Figs. 18 
to 20, illustrating the degree of agreement between simulated 
and observed peculiar velocity fields and between the real- 
space density fields reconstructed from them. For the simula- 
tions, the reconstructed density fields can also be compared 
with the real density fields which are plotted in Fig. 10. 
Overall the level of agreement is impressive. Systematic dif- 



ferences are visible in certain areas and presumably reflect 
noise elfects together with systematics arising because the 
smoothed IRAS 1.2 Jy density field which constrained our 
initial conditions does not perfectly trace the true mass den- 
sity field. Differences may also arise because the observed 
peculiar velocities are generated in part by the mass distri- 
bution outside the region where we constrained our initial 
conditions. 

Similar mock catalogues provided the basis for a re- 
vised likelihood analysis of peculiar velocities by Silberman 
et al. (2001). The proper incorporation of nonlinear effects 
revealed a systematic overestimate of Sim in earlier linear 
analyses. It then allowed the development of unbiased non- 
linear methods, which brought the best estimates from the 
Mark III and SFI data to Qm = 0.35 ±0.1. 

We publically release the mock catalogues for Mark 
III. They are available upon request from AE or AD (el- 
dar@phys.huji.ac.il). We can also provide similar mock cat- 
alogues tailor-made for other peculiar-velocity catalogues 
upon request. 



9 SUMMARY AND CONCLUSIONS 

The goal of the current study has been to carry out physi- 
cally consistent simulations of the formation and evolution 
of the local galaxy population within the currently dominant 
CDM structure formation paradigm. We have attempted 
to reproduce not only the statistics of galaxy properties 
(luminosity functions, colour and morphology distributions, 
TuUy-Fisher relations...) and of galaxy clustering (correla- 
tion functions, cluster luminosities...) but also the actual 
spatial distribution of clusters and superclusters within 80 
hr'^ Mpc of the Milky Way. This is an ambitious undertak- 
ing and clearly we have been more successful in some aspects 
of it than in others. 

On small scales the statistics of the dark matter dis- 
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Figure 17. Distribution of the mock apparent magnitudes and 
redshifts compared to the real data of the Mark III catalogue. 
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Figure 19. Same plot as Fig. 18, but for the mock rCDM cat- 
alogue. The density field shown has been reconstructed from the 
mock catalogue. 
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Figure 18. The IRAS 1.2 Jy map of radial peculiar velocities 
in a slice about the Supergalactic plane. The contours show the 
density field reconstructed from the IRAS data and smoothed 
over 12 h'^ Mpc. 
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Figure 20. Same plot as Fig. 19, but for the mock ACDM cata- 
logue. 



tributions in our ACDM and rCDM simulations are a very 
good match to those expected. On scales beyond a Mpc or 
two, discrepancies are apparent which reflect both the rel- 
atively small volume simulated and the constraints we im- 
posed on the initial large-scale fluctuations. For both simu- 
lations the morphology of the dark matter structure on large 
scales is a good match to that of the observed 1.2 Jy IRAS 
sample which we used as a constraint. In addition, the most 
massive clusters agree quite well with observed clusters both 
in position and in mass. 



The distributions of properties of individual galaxies are 
in fair agreement with those observed. The TuUy-Fisher re- 
lations we flnd for our simulated spiral galaxies are, by con- 
struction, in excellent agreement with observation. Our dis- 
tributions of morphology and colour also resemble those ob- 
served and show the correct dependence on galaxy environ- 
ment. Our luminosity functions, on the other hand, show sig- 
nificant discrepancies with those emerging from the 2dF and 
SDSS surveys. This problem is substantially reduced from 
K99, the paper on which our modelling is based, because of 
a number of minor improvements we have introduced, but 
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it is still worse than found by several other groups carrying 
out similar modelling (e.g. Cole et al. 2000; Somervillc ct al. 
2001, SOO). We conclude that this disagreement probably 
reflects an inadequacy in our galaxy modelling, rather than 
any fundamental problem with the paradigm. 

The statistics of gala^xy clustering in our ACDM simu- 
lation is in good agreement with that seen both in the far- 
infrared selected PSCz survey and in the optically selected 
UZC survey, at least as inferred from the autocorrelation 
functions of the galaxy populations. For our tCDM simula- 
tion the agreement with the PSCz survey is still quite good, 
but the amplitude of optical galaxy correlations is about 
40% lower than in the UZC survey. 

This latter problem may be related to the fact that 
the galaxy populations of our simulated galaxy clusters arc 
in poor agreement with observation. Our brightest cluster 
galaxies are mostly about a magnitude too bright, while 
other bright cluster galaxies are about a magnitude too faint. 
Overall our clusters have mass-to-light ratios at least a fac- 
tor of 2 larger than observed. This difficulty is apparently 
a consequence of the details of our galaxy formation mod- 
elling. Diaferio et al. (2001), using the simulations of K99, 
find a similar problem with the brightest cluster member 
but find their other cluster galaxies to be slightly too bright, 
while SOO find cluster luminosity functions which agree well 
in shape with observation, but mass-to-light ratios which are 
somewhat high. SOO show the improvement in luminosity 
function shape and the reduction in central galaxy luminos- 
ity to be a consequence of the greatly improved resolution of 
their simulation. This allows galaxy merging to be followed 
explicitly rather than inserted "by hand" using a dynamical 
friction model. 

Point-by-point comparison of the smoothed density 
fields of mass and of optical- and FIR-selected galaxies in 
our two simulations, as well as comparison with the observed 
density field of PSCz galaxies, illustrates both the nature of 
the biases which relate these various components and the 
degree to which our techniques have been successful in re- 
producing the observed spatial distribution of star-forming 
galaxies. Our simulations produce a remarkably tight rela- 
tion between the densities of optically selected galaxies and 
of dark matter. Biases are present but are relatively weak, 
too weak in the tCDM case to be consistent with obser- 
vation, as already noted above. The fact that stax-forming 
galaxies avoid rich groups and clusters results both in more 
scatter and in a strong bias at high density in the relation be- 
tween FIR galEixy density and either mass or optical galaxy 
density. Comparing density fields between simulations or 
with the observations shows the strong relations introduced 
by our initial condition constraints as well as biases caused 
by differing amounts of dynamical evolution. The excellent 
agreement between our simulated FIR galaxy distributions 
and the PSCz data is best illustrated by the cross-correlation 
between the two which is within 20% of the autocorrelation 
of the PSCz data on scales above 5 h^^ Mpc. 

Finally our mock catalogues of peculiar velocity data 
agree well with the real Mark III data. The reconstructed 
density fields from the mock catalogues show a close corre- 
spondancc to each other, to the true density fields in the 
simulations, and to the reconstructed density field from the 
real data. This is both a reassuring demonstration that our 
techniques have achieved their primary goal, and a sign- 



post to the way our simulations may be used to calibrate 
quantitative scientific analysis of the local galaxy distri- 
bution. We will present more detailed applications in fu- 
ture work, and we release a variety of galaxy, halo and 
dark matter catalogues in order that others can use them 
also. The data are available at the URL: http://www.mpa- 
garching.mpg.de/NumCos/CR/index.html 
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